Boundary effects in reaction-diffusion processes 
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The effects of a boundary on reaction systems are examined in the framework of the general single- 
species reaction/coalescence process. The boundary naturally represents the reactants' container, 
but is applicable to exciton dynamics in a doped TMMC crystal. We show that a density excess, 
which extends into the system diffusively from the boundary, is formed in two dimensions and below. 
This implies a surprising result for the magnetisation near a fixed spin in the coarsening of the one- 
dimensional critical Ising model. The universal, dimensionally-dependent functional forms of this 
density excess are given by an exact solution and the field-theoretic renormalisation group. 
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Reaction systems have been given a good deal of at- 
tention in the literature over the last few decades. They 
are widespread in nature and methods developed in their 
analysis have an applicability that extends well beyond 
conventional chemical systems. Moreover, they provide 
excellent examples of dynamical, many-body statistical 
processes and can exhibit a variety of interesting effects 
such as spontaneous symmetry breaking and pattern 
formation Q. 

In this work, we show that the imposition of impen- 
etrable boundaries on a reaction system leads to a non- 
trivial spatial variation in the reactant density. Even 
though such a boundary would naturally represent the 
vessel holding the reactants, until now no such studies 
have appeared in the literature. In systems with many 
interacting degrees of freedom, boundaries often give rise 
to surface effects that penetrate far into the bulk. Our 
results demonstrate that in reaction systems such long- 
range effects are indeed present and have a high degree 
of universality. 

We choose the class ]|,|| of the general single-species 
reaction processes as the starting point for the study of 
these boundary effects in reaction systems, as it will pro- 
vide a basis for the analysis of more complex systems. 
This universality class, comprising the annihilating ran- 
dom walk A+A — > O, coalescing random walk A+A — ► A 
and any combination thereof, is fundamental in the theo- 
retical study of reaction systems and covers a broad range 
of physical phenomena. For example, the coalescence 
process is seen in the dynamics of excitonic annihilation 
reactions in the TMMC crystal [§§]. The predicted de- 
cay exponent from theory is in agreement with exper- 
iment for over five orders of magnitude. It should be 
noted that the boundary effects we introduce here are 
exhibited in the TMMC crystal in the presence of Mg 2+ 
or Cd 2+ doping. These defect ions act as perfect reflec- 
tors for the annihilating excitons j7j. 

A mapping also exists from the annihilating random 
walk to the domain coarsening dynamics in the critical 
one-dimensional Ising magnet |^ [lC[] . In this mapping 
the A particles represent domain walls with an annihila- 



tion of two A particles corresponding to a domain shrink- 
ing to zero size in a background of the opposite phase. 
An impenetrable boundary in the reaction system cor- 
responds to a fixed boundary spin in the Ising magnet. 
Our analysis will show that domain walls are more likely 
to be found near the fixed spin than far away. This gives 
the counter-intuitive result that the absolute value of the 
dynamic, coarse-grained magnetisation is actually lower 
near a fixed spin. 

In the interests of notational simplicity, we provide 
here the analysis specifically for the process A + A — > O. 
However, all densities given can be converted to the result 
for the coalescence process A + A — > A by a simple factor 
of two. The annihilating random walk has been studied 
extensively in homogeneous unbounded systems, cither 
of infinite extent or with periodic boundary conditions. 
Throughout the work, this will be referred to as the bulk 
case. It is well known ||| 11 1^] that the variation of the 
density g with time t differs from the mean-field predic- 
tion q oc 1/f for dimensions d < 2. In fact the actual 
density decays are g oc t~ x logt id — 2), and g oc t~ d / 2 
(d < 2) with a universal amplitude. In all cases, it must 
be stressed that the density remains uniform throughout 
the system. 

In the following, we will first introduce the boundary 
into the reaction process and specify the model to be 
studied. The mean-field approximation will be shown 
to predict a homogeneous density unchanged from the 
bulk case. However, an argument will be presented to 
show that the mean-field prediction breaks down in low 
dimensions. The central result of this work is that a 
fluctuation-induced excess density develops at the bound- 
ary and extends into the system diffusively. We outline 
the field-theoretic renormalisation group (RG) descrip- 
tion which we use to identify the universal quantities of 
this excess. These calculations were performed in real 
space to one-loop order and show that in two dimension 
and below, the density excess gE has the following form 
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Here, D is the diffusion constant of the reactants, z is 
the normal distance from the boundary and fd is the di- 
mensionally dependent scaling function. We were able 
to find the late-time scaling functions ji and /i exactly 
in both two and one dimensions (given in Eqs. (|^) and 
( pd| ) respectively). The former we derive from the RG 
improved field-theoretic calculation and the latter from 
an exact solution. Finally, in the context of surface criti- 
cal phenomena, the behaviour in the related ballistic an- 
nihilation reaction system will be briefly examined and 
compared with the diffusive case. 

We now introduce the model. The system is defined on 
a hypercubic, d-dimensional lattice with a lattice spac- 
ing of unity. The lattice, is infinite in d — 1 transverse 
dimensions and semi-infinite (sites 1,2, •••,00) in what 
will be called the z direction. At time t = the lattice 
is filled with an initial density Qq of identical particles 
that perform two types of dynamics: diffusion and mu- 
tual annihilation. The diffusion is represented by each 
particle hopping at a rate D to any neighbouring site at 
random. A hop from site z — 1 towards the boundary is 
disallowed. The diffusion is independent for each parti- 
cle, and hence multiple occupancy of each site is possible, 
leading to bosonic particle statistics. However, if there 
are n > 2 particles on a particular site, a reaction can 
occur there with a rate Xn(n — 1) to reduce n by 2. 

The above dynamics can be approximated by a mean- 
field description. This involves ignoring all possible cor- 
relations by considering a self-consistent equation involv- 
ing just the average density ~g. In the continuum limit the 
boundary is enforced by a zero-current restriction, thus 



d t g = DV 2 g - 2Xg 2 with d z g\ z=a = 0. 
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The boundary restriction is compatible with the bulk so- 
lution ~g = go/(l + go2Xt). Hence, in the absence of strong 
fluctuations, the density is uniform throughout the sys- 
tem. 

However, correlations must be properly accounted for 
in low dimensions. Consider the dynamics of the model, 
in one dimension, up to a time t and far from the wall. 
Because random walks in one-dimension are recurrent, 
most particles within a diffusion length lb ~ \/2Dt in the 
bulk, will have interacted and annihilated. This leads to 



a density in the bulk of the system of gb — 1^ — Cb/y/t. 
However, close to the wall the diffusion length is smaller, 
g w ~ l^ 1 ~ c w j\ft since c w > Cb, leading to a density 
excess near the boundary. 

Nevertheless, this argument is rather crude and a 
method for systematically including fluctuations is re- 
quired. The RG has provided such a method for calcu- 
lating bulk quantities P, pd|Jl3| - [l5| ] , with the advantage of 
clearly identifying universal properties. We now present 
an overview of the generalisation to a system with a 
boundary: a case technically more complicated due to 
the lack of translational invariance. Details of the calcu- 
lation will be provided elsewhere [pi . 



The field-theoretic description is obtained by first writ- 
ing a master equation. This describes the flow of prob- 
ability between microstates of the system and is con- 
veniently written in second-quantised form dt\P >= 
—H\P>. The vector |P> is the probability-state vector 
written in a Fock space and acted upon by the evolution 
operator TL 
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where a\a are the usual bosonic operators. The sum i 
is over all lattice sites, and the sum j is over all of site 
i's neighbours, with the condition that both sums are 
restricted to the half-space. 

The algebraic description 7i{a, (v) is mapped onto a 
continuum path-integral for the action S(<f), <fi) where the 
complex fields 4> and <j> are analogous to a and a' —1. The 
action S = Sd + Sr + S eo comprises diffusion So, reac- 
tion Sr and initial-condition S go components. The diffu- 
sive part provides the propagator for the theory which is 
Gaussian for the transverse dimensions and has the the 
following "mirrored" form for the z-dimension 



Gz(zf, Zi,t) = G(zf - Zi,t) + G(zf + Zi,t) 



(4) 



where G(z, i) is a Gaussian with a standard deviation of 
2Dt. The other components in the action are 



S n 2 A / dt f d d r<j)(j) 2 + X f dt [ d d r(f/ > 2 , (5) 

Jz>0 JO Jz>Q 

t 



S eo =-Qo dt d d r<j>S(t). 

JO J 2>0 

The upper critical dimension of the theory is d c = 2 and 
observables were rendered finite by dimensional regurali- 
sation in d = 2 — e. The propagator is not dressed by the 
interactions (|^), implying the boundary remains effec- 
tively reflecting on all scales. In the language of surface 
critical phenomena, this corresponds to the special tran- 
sition fill persisting at all orders. This is different from 
the behaviour frequently seen in equilibrium surface crit- 
ical phenomena |l8f| and in related non-equilibrium sys- 
tems [fl9f . In fact only the reaction rate A is renormalised, 
with a fixed point structure identical to the bulk case jll] . 
This is understandable as physically the renormalisation 
of A is connected to the fact that random walks in two 
dimensions and below are recurrent: a feature unaffected 
by the presence of a boundary. 

To get non-trivial, z-dependent results it is clear (from 
the lack of an excess in the mean-field equation) that 
the RG improved perturbation expansion must be taken 
to at least one-loop order. Writing the density as an 
expansion in e = 2 — d and splitting the contribution into 
an excess qe{z,£) and a homogeneous, background bulk 
density £>s(i) 
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g(z,t) = g B (t) + g E (z,t), 



(6) 



the homogeneous bulk density for d < 2 is found to be 
1 
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This is exactly the result found in |Tl]] as expected. 

However, a fluctuation-induced density excess is also 
found, representing the new result from this calculation 
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The function EfJ has the asymptotics ~ exp 
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and is plotted in Fig. 1. A few things should be noted 
about the form of g E . First, this excess is not localised at 
the boundary but extends into the system diffusively, by 
virtue of its functional dependence of z 2 /Dt. Also, the 
excess shares the same universality as the bulk density 
in that it is independent of the reaction rate A and the 
initial density g$. Finally, for d < 2 the amplitude of 
the excess decays with the same exponent as the bulk. 
Hence, the ratio of the boundary to the bulk density 
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is a constant, universal quantity independent of all sys- 
tem parameters except the dimension. 

The behaviour in two dimensions provides for an in- 
teresting result: the one-loop calculation Eq. (||) is the 
exact late-time density excess. It is independent from 
the renormalised reaction rate and therefore represents 
the universal leading order, with higher-loop corrections 
decaying as (ilogi) -1 . The excess given in Eq. (||) gives 
surprisingly accurate results even for short times, as can 
be seen in Fig. 1 where the result for the density excess 
is compared with simulations. 

The RG has provided information about the general 
behaviour of the density excess as a function of dimen- 
sion. The universal quantities have been identified and 
the excess density correctly predicted for late times in 
d = 2. Unfortunately, the e expansion gives disappoint- 
ing results for the amplitude in d = 1 . Motivated by the 
mapping of the reaction dynamics onto the Ising magnet 
and the related excitonic coalescence process, we now 
provide the exact solution in one dimension. We gen- 
eralise the bulk model described in to include an 
impenetrable boundary. The on-site reaction rate is cho- 
sen to be infinite and hence a site in the system can be 
occupied by at most one particle. The dynamical rules 
can be written in terms of the possible evolutions of a 
pair of neighbouring sites. Denoting a particle on site k 
by Ak and an empty site by Ok, the allowed changes are 




1 
scaled distance from the boundary z/(2Dt) m 

FIG. 1. The density excess in one and two dimensions, /i 
and f%. The RG result fa is compared with simulations at 
time t — 10 (o) and t = 80 (*). The exact result /i calculated 
for infinite reaction rate is compared with a simulation for fi- 
nite reaction rate A = 1/2 (+) demonstrating the universality. 



O z A z+ i <-> A z O z +i with a rate D 
A z A z+ i — > O z O z +i with a rate 2D 
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where the site label z — 1 , 2 • ■ • is restricted to positive 
integers. The master equation can be written in the lan- 
guage of spin-half operators, and we were able to solve 
for the density using a similarity transformation found 
in |H| (details of the calculation will be presented else- 
where p6|). The density can be written as a sum of 
Bessel functions, of which the continuum limit gives the 
following density excess 
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h(e) = V2erfc exp 



erfc (£) . 



This function has the asymptotics /i ~ exp~£ I 2 /£ and 
is plotted in Fig. 1. Again, the bulk component of the 
density g E = l/V87rDt was found to be identical to the 
infinite case p2[ . As expected, the form of the solution is 
in agreement with that predicted from the RG treatment. 
The universal ratio of the density near the boundary to 
the bulk density given in Eq. (^) is found to be exactly 
\/2. It should be noted that the one-loop RG prediction 
for the universal ratio (given in Eq. (JsJ) with e = 1) is 
numerically ~ 1.47 for one dimension, and therefore gives 
quite a fair indication of the exact value. This should 
be compared with the poor one-loop predictions for the 
amplitudes of the density itself in one dimension. 

The universality (independence from the reaction rate 
A) predicted by the RG treatment, can be seen by com- 
paring this exact result for A = oo with data from a 
simulation of a system with finite reaction rate (Fig. 1). 
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The result, Eq. (|ll"|), gives the time-dependent prob- 
ability density of domain walls in the coarsening, one- 
dimensional critical Ising model near a very strong mag- 
netic field, i.e. a fixed spin. The magnitude of the coarse- 
grained magnetisation is a function of the local density 
of domain walls - the fewer the domain walls the higher 
the magnetisation. Hence, the result implies that the ab- 
solute value of the magnetisation measured near a fixed 
spin is lower than in the bulk of the system. It would 
be very interesting to see if this dynamic effect is seen in 
other magnetic systems or in higher dimensions. 

In summary, results have been presented from an anal- 
ysis of a reaction-diffusion process near an impenetrable 
boundary. The mean-field equation was shown to pre- 
dict a flat density profile. However, it was demonstrated 
that in two dimensions and below, a density excess de- 
velops at the boundary. In one dimension it was found 
that the density excess is as significant a contribution as 
the bulk density with both decaying as ~ t^ 1 / 2 . In two 
dimensions the excess was found to be marginally sub- 
dominant, decaying as t . Both these density excesses 
shares the same universality as the bulk density, in that 
they are independent of the reaction rate and the initial 
density. Moreover, a higher degree of universality was 
found in the ratio between the boundary and bulk den- 
sities: a quantity depending only on the dimension of 
space. The functional forms for the density excess were 
obtained for both one and two dimensions and a link was 
made to the coarsening dynamics of the one-dimensional 
Ising model. Most importantly, it was shown that the 
excess is not localised at the boundary but extends dif- 
fusively throughout the bulk of the system. It should 
be noted that above the critical dimensions there is a 
sub-dominant density excess. However, these effects are 
transient and quickly decay to yield the mean-field result, 
as predicted. 

A quantity much studied recently and related to do- 
main coarsening dynamics is the persistence exponent. 
This exponent describes the time dependence of the dis- 
tribution of sites not yet visited by a domain wall. In the 
homogeneous bulk case, persistence has been examined 
by the RG J22|,^3| and an exact result found in one di- 
mension p4] . In light of result ( |TT| ) and the fact that the 
excess penetrates into the system, it would be worth ex- 
amining the behaviour of the persistence in the presence 
of a fixed spin. 

Finally, it is interesting to consider the difference be- 
tween the behaviour of diffusing and gas-phase ballistic 
reaction kinetics near a boundary. A simple realisation 
of the ballistic A+A-^O system shows the same bulk 
density-decay exponent in one dimension as the diffusive 
case. However, contrary to the diffusive case, it can be 
shown that for ballistic reactions with an (elastic) im- 
penetrable boundary there is a lower density of reactants 
near the wall fig ]. This is understandable, because in 
late times the remaining particles tend to congregate in 



groups moving in the same direction. When such a group 
hits the elastic wall it mostly annihilates within itself, 
leaving few particles to return. In the context of surface 
critical phenomena this corresponds to the impenetrable 
boundary behaving as an effectively absorbing boundary 
on long time scales, in contrast to the diffusive case. 
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